#define MAX_SIZE 1111
#define MAX_THREADS 8
*****************************************************************************
*									    *
*	1D (slow) Fourier Transform					    *
*									    *
*****************************************************************************
      subroutine dft1du( plusminus, n, a, inc)
      double precision a(*)
      double precision w(MAX_SIZE)
      integer n, plusminus, inc
      
      double precision u, v, tpi
      TPI = 2 * 3.14159265358979323846

      if( inc .ne . 1) goto 11
      if ( plusminus .eq. -1) then
	l = (n+2)/2
c$doacross local(k,u,v), share(w)
	do k = 1, l
	    u = 0.0
	    v = 0.0
	    do i = 1, n
		u = u + a(i) * cos((k-1) * (i-1) * tpi / float(n))
		v = v - a(i) * sin((k-1) * (i-1) * tpi / float(n))
	    end do
	    w((2*k-2)+1) = u
	    w((2*k-1)+1) = v
	end do
c$doacross
	do i = 1, 2*l
	    a(i) = w(i)
	end do
      elseif (plusminus .eq. 1) then
	print *,'NOT Implemented'
      else
	print *, 'Error first argument of DFT1DU should -1 or +1'
	stop
 11	print *, 'This Simple version support only INC = 1'
	stop
      endif

      return
      end


      subroutine zft1d( plusminus, n, a, inc)
      implicit double precision ( a-h, o-z)
      double complex a(*)
      double complex w(MAX_SIZE)
      integer n, plusminus
      
      double complex cp
      double precision re, im, tpi

      tpi = 2.*3.14159265358979323846

	do i = 1, n
	    w(i) = 0.0
	    do k = 1, n
		im = mod(((i-1)*(k-1)), n)
		im = dfloat( plusminus*im)*tpi / dfloat(n) 
		re = cos( im)
		im = sin( im)
		cp = dcmplx( re, im)
		w(i) = w(i) + cp*a(inc*(k-1)+1)
	    end do
	end do
	do i = 1, n
	    a(inc*(i-1)+1) = w(i)
	end do

      return
      end



*****************************************************************************
*									    *
*	2D (slow) Fourier Transform					    *
*									    *
*****************************************************************************
      subroutine zft2d( plusminus, n1, n2, w, ldw)
      double complex w(ldw,n2)
      double complex save( MAX_SIZE)
      integer n1, n2, lda, ldw, plusminus

      call zfft1di( n1, save)
c$doacross local(j)
      do j = 1, n2
	call zfft1d( plusminus, n1, w(1,j), 1, save)
      end do

      call zfft1di( n2, save)
c$doacross local(i,j,k,re,im,cp)
      do i = 1, n1
	call zfft1d( plusminus, n2, w(i,1), ldw, save)
      end do

      return
      end

      subroutine dft2du( plusminus, n1, n2, w, ldw)
      Double Precision w(ldw,n2)
      Double Precision save( 2*MAX_SIZE), cp(2*MAX_SIZE,MAX_THREADS)
      integer n1, n2, lda, ldw, plusminus, l
 
      l = ((n1+2)/2)*2
      if (ldw .lt. l) then
	PRINT *, 'ERROR Leading Dimension is not sufficient',short(l),short(ldw)
      endif
      if (plusminus .eq. -1) then
	call dfft1di( n1, save)
c$doacross local(j)
	do j = 1, n2
	    call dfft1du( plusminus, n1, w(1,j), 1, save)
	end do

	call zfft1di( n2, save)
c$doacross local(i,j,my), share(cp)
	do i = 1, l, 2
	    my = 1
c$	    my = mp_my_threadnum()+1
	    do j = 1, n2
		cp(2*j-1,my) = w(i+0,j)
		cp(2*j-0,my) = w(i+1,j)
	    end do
	    call zfft1d( plusminus, n2, cp(1,my), 1, save)
	    do j = 1, n2
		w(i+0,j) = cp(2*j-1,my) 
		w(i+1,j) = cp(2*j-0,my)
	    end do
	end do
      else
	call zfft1di( n2, save)
c$doacross local(i,j,my)
	do i = 1, l, 2
	    my = 1
c$	    my = mp_my_threadnum()+1
	    do j = 1, n2
		cp(2*j-1,my) = w(i+0,j)
		cp(2*j-0,my) = w(i+1,j)
	    end do
	    call zfft1d( plusminus, n2, cp(1,my), 1, save)
	    do j = 1, n2
		w(i+0,j) = cp(2*j-1,my) 
		w(i+1,j) = cp(2*j-0,my)
	    end do
	end do

	call dfft1di( n1, save)
c$doacross local(j)
	do j = 1, n2
	    call dfft1du( plusminus, n1, w(1,j), 1, save)
	end do
      end if
      return
      end
      

*****************************************************************************
*
*	3D (slow) Fourier Transform
*
*****************************************************************************

      subroutine zft3d(plusminus,n1,n2,n3,w,ld1,ld2)
      double complex w(ld1,ld2,n3)
      double complex save( MAX_SIZE)
      integer n1,n2,n3,ld1,ld2,plusminus
c
c   transform along X  first ...
c
      call zfft1di( n1, save)
c$doacross local(i,j,k)
      do k = 1, n3
       do j = 1, n2
	    call zfft1d( plusminus, n1, w(1,j,k), 1, save)
       end do
      end do
c
c   transform along Y then ...
c
      call zfft1di( n2, save)
c$doacross local(i,j,k)
      do k = 1,n3
       do i = 1,n1
	    call zfft1d( plusminus, n2, w(i,1,k), ld1, save)
       end do
      end do
c
c   transform along Z finally ...
c
      call zfft1di( n3, save)
c$doacross local(i,j,k)
      do i = 1, n1
       do j = 1, n2
	    call zfft1d( plusminus, n3, w(i,j,1), ld1*ld2, save)
       end do
      end do

      return
      end

      subroutine dft3du(plusminus,n1,n2,n3,w,ld1,ld2)
      Double Precision w(ld1,ld2,n3)
      Double Precision save(2*MAX_SIZE), cp(MAX_SIZE,MAX_THREADS)
      integer n1,n2,n3,ld1,ld2,plusminus, l

      if( plusminus .ne. -1) then
	print *, 'Option non implemented YET'
      end if
      l = ((n1+2)/2)*2
      if( ld1 .lt. l ) then
	print *, 'Leading dimension not sufficient'
      endif
c
c   trandform along X  first ...
c
      call dfft1di( n1, save)
c$doacross local(i,j,k)
      do k = 1, n3
       do j = 1, n2
	    call dfft1du( plusminus, n1, w(1,j,k), 1, save)
       end do
      end do
c
c   trandform along Y then ...
c
      call zfft1di( n2, save)
c$doacross local(i,j,k,my), share(cp)
      do k = 1,n3
	my = 1
c$	my = mp_my_threadnum()+1
	do i = 1, l, 2
	    do j = 1, n2
		cp(2*j-1,my) = w(i+0,j,k)
		cp(2*j-0,my) = w(i+1,j,k)
	    end do
	    call zfft1d( plusminus, n2, cp(1,my), 1, save)
	    do j = 1, n2
		w(i+0,j,k) = cp(2*j-1,my) 
		w(i+1,j,k) = cp(2*j-0,my)
	    end do
	end do
      end do
c
c   trandform along Z finally ...
c
      call zfft1di( n3, save)
c$doacross local(i,j,k,my)
      do j = 1,n2
	my = 1
c$	my = mp_my_threadnum()+1
	do i = 1, l, 2
	    do k = 1, n3
		cp(2*k-1,my) = w(i+0,j,k)
		cp(2*k-0,my) = w(i+1,j,k)
	    end do
	    call zfft1d( plusminus, n3, cp(1,my), 1, save)
	    do k = 1, n3
		w(i+0,j,k) = cp(2*k-1,my)
		w(i+1,j,k) = cp(2*k-0,my)
	    end do
	end do
      end do

      return
      end

*****************************************************************************
*									    *
*	1D (slow) Fourier Transform					    *
*									    *
*****************************************************************************
      subroutine sft1du( plusminus, n, a, inc)
      real a(*)
      real w(MAX_SIZE)
      integer n, plusminus, inc
      
      real u, v, tpi
      TPI = 2 * 3.14159265358979323846

      if( inc .ne . 1) goto 11
      if ( plusminus .eq. -1) then
	l = (n+2)/2
c$doacross local(k,u,v), share(w)
	do k = 1, l
	    u = 0.0
	    v = 0.0
	    do i = 1, n
		u = u + a(i) * cos((k-1) * (i-1) * tpi / float(n))
		v = v - a(i) * sin((k-1) * (i-1) * tpi / float(n))
	    end do
	    w((2*k-2)+1) = u
	    w((2*k-1)+1) = v
	end do
c$doacross
	do i = 1, 2*l
	    a(i) = w(i)
	end do
      elseif (plusminus .eq. 1) then
 11	print *,'NOT Implemented'
      endif

      return
      end


      subroutine cft1d( plusminus, n, a, inc)
      implicit real ( a-h, o-z)
      complex a(*)
      complex w(MAX_SIZE)
      integer n, plusminus
      
      complex cp
      real re, im, tpi

      tpi = 2.*3.14159265358979323846

	do i = 1, n
	    w(i) = 0.0
	    do k = 1, n
		im = mod(((i-1)*(k-1)), n)
		im = float( plusminus*im)*tpi / float(n) 
		re = cos( im)
		im = sin( im)
		cp = cmplx( re, im)
		w(i) = w(i) + cp*a(inc*(k-1)+1)
	    end do
	end do
	do i = 1, n
	    a(inc*(i-1)+1) = w(i)
	end do

      return
      end



*****************************************************************************
*									    *
*	2D (slow) Fourier Transform					    *
*									    *
*****************************************************************************
      subroutine cft2d( plusminus, n1, n2, w, ldw)
      complex w(ldw,n2)
      complex save( MAX_SIZE)
      integer n1, n2, lda, ldw, plusminus

      call cfft1di( n1, save)
c$doacross local(j)
      do j = 1, n2
	call cfft1d( plusminus, n1, w(1,j), 1, save)
      end do

      call cfft1di( n2, save)
c$doacross local(i,j,k,re,im,cp)
      do i = 1, n1
	call cfft1d( plusminus, n2, w(i,1), ldw, save)
      end do

      return
      end

      subroutine sft2du( plusminus, n1, n2, w, ldw)
      real w(ldw,n2)
      real save( 2*MAX_SIZE), cp(2*MAX_SIZE,MAX_THREADS)
      integer n1, n2, lda, ldw, plusminus, l
 
      l = ((n1+2)/2)*2
      if (ldw .lt. l) then
	PRINT *, 'ERROR Leading Dimension is not sufficient',short(l),short(ldw)
      endif
      if (plusminus .eq. -1) then
	call sfft1di( n1, save)
c$doacross local(j)
	do j = 1, n2
	    call sfft1du( plusminus, n1, w(1,j), 1, save)
	end do

	call cfft1di( n2, save)
c$doacross local(i,j,my), share(cp)
	do i = 1, l, 2
	    my = 1
c$	    my = mp_my_threadnum()+1
	    do j = 1, n2
		cp(2*j-1,my) = w(i+0,j)
		cp(2*j-0,my) = w(i+1,j)
	    end do
	    call cfft1d( plusminus, n2, cp(1,my), 1, save)
	    do j = 1, n2
		w(i+0,j) = cp(2*j-1,my) 
		w(i+1,j) = cp(2*j-0,my)
	    end do
	end do
      else
	call cfft1di( n2, save)
c$doacross local(i,j,my)
	do i = 1, l, 2
	    my = 1
c$	    my = mp_my_threadnum()+1
	    do j = 1, n2
		cp(2*j-1,my) = w(i+0,j)
		cp(2*j-0,my) = w(i+1,j)
	    end do
	    call cfft1d( plusminus, n2, cp(1,my), 1, save)
	    do j = 1, n2
		w(i+0,j) = cp(2*j-1,my) 
		w(i+1,j) = cp(2*j-0,my)
	    end do
	end do

	call sfft1di( n1, save)
c$doacross local(j)
	do j = 1, n2
	    call sfft1du( plusminus, n1, w(1,j), 1, save)
	end do
      end if
      return
      end
      

*****************************************************************************
*
*	3D (slow) Fourier Transform
*
*****************************************************************************

      subroutine cft3d(plusminus,n1,n2,n3,w,ld1,ld2)
      complex w(ld1,ld2,n3)
      complex save( MAX_SIZE)
      integer n1,n2,n3,ld1,ld2,plusminus
c
c   transform along X  first ...
c
      call cfft1di( n1, save)
c$doacross local(i,j,k)
      do k = 1, n3
       do j = 1, n2
	    call cfft1d( plusminus, n1, w(1,j,k), 1, save)
       end do
      end do
c
c   transform along Y then ...
c
      call cfft1di( n2, save)
c$doacross local(i,j,k)
      do k = 1,n3
       do i = 1,n1
	    call cfft1d( plusminus, n2, w(i,1,k), ld1, save)
       end do
      end do
c
c   transform along Z finally ...
c
      call cfft1di( n3, save)
c$doacross local(i,j,k)
      do i = 1, n1
       do j = 1, n2
	    call cfft1d( plusminus, n3, w(i,j,1), ld1*ld2, save)
       end do
      end do

      return
      end

      subroutine sft3du(plusminus,n1,n2,n3,w,ld1,ld2)
      real w(ld1,ld2,n3)
      real save(2*MAX_SIZE), cp(MAX_SIZE,MAX_THREADS)
      integer n1,n2,n3,ld1,ld2,plusminus, l

      if( plusminus .ne. -1) then
	print *, 'Option non implemented YET'
      end if
      l = ((n1+2)/2)*2
      if( ld1 .lt. l ) then
	print *, 'Leading dimension not sufficient'
      endif
c
c   transform along X  first ...
c
      call sfft1di( n1, save)
c$doacross local(i,j,k)
      do k = 1, n3
       do j = 1, n2
	    call sfft1du( plusminus, n1, w(1,j,k), 1, save)
       end do
      end do
c
c   transform along Y then ...
c
      call cfft1di( n2, save)
c$doacross local(i,j,k,my), share(cp)
      do k = 1,n3
	my = 1
c$	my = mp_my_threadnum()+1
	do i = 1, l, 2
	    do j = 1, n2
		cp(2*j-1,my) = w(i+0,j,k)
		cp(2*j-0,my) = w(i+1,j,k)
	    end do
	    call cfft1d( plusminus, n2, cp(1,my), 1, save)
	    do j = 1, n2
		w(i+0,j,k) = cp(2*j-1,my) 
		w(i+1,j,k) = cp(2*j-0,my)
	    end do
	end do
      end do
c
c   transform along Z finally ...
c
      call cfft1di( n3, save)
c$doacross local(i,j,k,my)
      do j = 1,n2
	my = 1
c$	my = mp_my_threadnum()+1
	do i = 1, l, 2
	    do k = 1, n3
		cp(2*k-1,my) = w(i+0,j,k)
		cp(2*k-0,my) = w(i+1,j,k)
	    end do
	    call cfft1d( plusminus, n3, cp(1,my), 1, save)
	    do k = 1, n3
		w(i+0,j,k) = cp(2*k-1,my)
		w(i+1,j,k) = cp(2*k-0,my)
	    end do
	end do
      end do

      return
      end

